function handles = analyzeInputMapArray(handles)
% analyzeInputMapArray     Main engine for analyzing set of traces.
%
%
% gs april 2005 - Modified for mapAnalysis2p0 from old analyzeArray
% ------------------------------------------------------------------------------

% ANALYSIS METHOD
analysisMethod = get(handles.methodInputMap, 'Value');
% method 1 -- old way: baseline-subtract using MEAN; divide baseline into baseline/response
% method 2 -- new way: baseline-subtract using MEDIAN; simply calc SD of baselines

% DATA etc
dataArray = handles.data.map.mapActive.bsArray; % dataArray
% should modify this to make it conditional ...

[rows, numTraces] = size(dataArray); % must be 2D data array; col num = trace num
sr = handles.data.acq.sampleRate;
mapPattern = handles.data.map.mapActive.mapPattern;
traceData = [];

% CELL parameters
if get(handles.chkWithRsRm, 'Value')
    if ~handles.data.analysis.currentClamp
        handles.data.analysis.skipCm = 0; % option here to skip Cm calc, e.g. if it's too slow
        handles = arrayCellParams(handles, dataArray);
    else
        disp('Cell parameter calculations not yet implemented for current clamp.');
    end
end

% TIME vector
totalTime = (rows-1)/sr; % in sec, not msec
xTimeAxis = linspace(0, totalTime, rows)';
handles.data.map.mapActive.xTimeAxis = xTimeAxis;
minX = 0; 
maxX = max(xTimeAxis);

% BASELINE settings, mean, sd, baseline-subtracted array
if analysisMethod == 1
    baselineStartIndex = round(handles.data.analysis.baselineStart*sr)+1;
    baselineEndIndex = round(handles.data.analysis.baselineEnd*sr)+1;
    baselineMeans = mean(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineMeans = baselineMeans';
    baselineSDs = std(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineSDs = baselineSDs';
    % % baseline-subtracted array:
    % baselineArray = repmat(baselineMeans, rows, 1);
    % bsArray = dataArray - baselineArray; 
    % handles.data.map.mapActive.bsArray = bsArray;
    bsArray = handles.data.map.mapActive.bsArray;
    
    % NOISE of baseline
    %       b0 = 0-50 msec, baseline for baseline window
    %       r0 = 50-100 msec, baseline window
    %       b1 = 50-100 msec, baseline for synaptic window
    %       r1 = 110-160 msec, synaptic window
    % % mean 
    % b0 = mean(state.uncaging.map.bsArray(1:500, :));
    % b1 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r0 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r1 = mean(state.uncaging.map.bsArray(1101:1601, :));
    % peak 
    b0 = min(handles.data.map.mapActive.bsArray(1:500, :));
    b1 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r0 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r1 = min(handles.data.map.mapActive.bsArray(1101:1601, :));
    % 
    handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);
    handles.data.map.mapActive.histoBaseSDs = std(handles.data.map.mapActive.ampHistoBaselineData);
    handles.data.map.mapActive.ampHistoSynapticData = -(r1-b1);
elseif analysisMethod == 2
    baselineStartIndex = round(handles.data.analysis.baselineStart*sr)+1;
    baselineEndIndex = round(handles.data.analysis.baselineEnd*sr)+1;
    baselineMedians = median(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineMedians = baselineMedians';
    
    baselineSDs = std(dataArray(baselineStartIndex:baselineEndIndex,:));
    handles.data.map.mapActive.baselineSDs = baselineSDs';
    
    % baseline-subtracted array:
    baselineArray = repmat(baselineMedians, rows, 1);
    bsArray = dataArray - baselineArray; 
    handles.data.map.mapActive.bsArray = bsArray;
    
    % NOISE of baseline
    %       b0 = 0-50 msec, baseline for baseline window
    %       r0 = 50-100 msec, baseline window
    %       b1 = 50-100 msec, baseline for synaptic window
    %       r1 = 110-160 msec, synaptic window
    % % mean 
    % b0 = mean(state.uncaging.map.bsArray(1:500, :));
    % b1 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r0 = mean(state.uncaging.map.bsArray(500:1000, :));
    % r1 = mean(state.uncaging.map.bsArray(1101:1601, :));
    % peak 
    b0 = min(handles.data.map.mapActive.bsArray(1:500, :));      % these are all vectors
    b1 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r0 = min(handles.data.map.mapActive.bsArray(500:1000, :));
    r1 = min(handles.data.map.mapActive.bsArray(1101:1601, :));
    % 
    handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);  % vector
    handles.data.map.mapActive.histoBaseSDs = std(handles.data.map.mapActive.ampHistoBaselineData); % single value (not a vector)
    
    % new way
    handles.data.map.mapActive.ampHistoBaselineData = -(r0-b0);
    
    handles.data.map.mapActive.histoBaseSDs = mean(std(dataArray(baselineStartIndex:baselineEndIndex,:)));
    handles.data.map.mapActive.ampHistoSynapticData = -(r1-b1);
end

% STIMULUS timing
stimOn = handles.data.analysis.stimOn;
stimOnInd = round(stimOn*sr+1);
handles.data.analysis.stimOnInd = stimOnInd;

% WINDOWS *******************************************************************
directWindowStart = handles.data.analysis.directWindowStart;
directStartInd = round(directWindowStart * sr + 1);
directWindow = handles.data.analysis.directWindow;
directEndInd = directStartInd + round(directWindow * sr);

synapticWindowStart = handles.data.analysis.synapticWindowStart;
synapticWindowStartInd = round(synapticWindowStart * sr + 1);
synapticWindow = handles.data.analysis.synapticWindow;
synEndInd = directStartInd + round(synapticWindow * sr);

fourthWindowStart = handles.data.analysis.fourthWindowStart;
fourthWindowStartInd = round(fourthWindowStart * sr + 1);
fourthWindow = handles.data.analysis.fourthWindow;
fourthWindowEndInd = fourthWindowStartInd + round(fourthWindow * sr);

synDuration = handles.data.analysis.synDuration;

% THRESHOLDS ****************************************************************
% NB: these are vectors
synLevel = handles.data.analysis.synThreshold; 
dirLevel = synLevel; 
if analysisMethod == 1 % method 1 -- based on baseline mean
    handles.data.map.mapActive.dirPosThresh = handles.data.map.mapActive.baselineMeans + dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.dirNegThresh = handles.data.map.mapActive.baselineMeans - dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synPosThresh = handles.data.map.mapActive.baselineMeans + synLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synNegThresh = handles.data.map.mapActive.baselineMeans - synLevel * handles.data.map.mapActive.histoBaseSDs;
elseif analysisMethod == 2 % method 1 -- based on baseline median
    handles.data.map.mapActive.dirPosThresh = dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.dirNegThresh = -dirLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synPosThresh = synLevel * handles.data.map.mapActive.histoBaseSDs;
    handles.data.map.mapActive.synNegThresh = -synLevel * handles.data.map.mapActive.histoBaseSDs;
end

handles.data.map.mapActive.synThreshpAmV = synLevel * handles.data.map.mapActive.histoBaseSDs;
handles.data.map.mapActive.baselineSD = handles.data.map.mapActive.histoBaseSDs;

% RESPONSES *****************************************************************
% strategy: measure response params over entire post-stim window (direct + syn)
% aftwerward, determine whether responses were direct or syn

% MEAN
genMeanAmp = mean(bsArray(stimOnInd:synEndInd,:));
handles.data.map.mapActive.genMeanAmp = genMeanAmp';

% MIN
[genMinAmp, genMinPkLatInd] = min(bsArray(stimOnInd:synEndInd,:));
genLatOfMin = (genMinPkLatInd-1)/sr;
handles.data.map.mapActive.genMinAmp = genMinAmp';
handles.data.map.mapActive.genLatOfMin = genLatOfMin';

% MAX
[genMaxAmp, genMaxPkLatInd] = max(bsArray(stimOnInd:synEndInd,:));
genLatOfMax = (genMaxPkLatInd-1)/sr;
handles.data.map.mapActive.genMaxAmp = genMaxAmp';
handles.data.map.mapActive.genLatOfMax = genLatOfMax';

% INTEGRAL
genIntegral = trapz(bsArray(stimOnInd:synEndInd,:))/sr;
handles.data.map.mapActive.genIntegral = genIntegral';

% 'FOURTH WINDOW'
handles.data.map.mapActive.genMeanAmp4th = mean(bsArray(fourthWindowStartInd:fourthWindowEndInd,:));
handles.data.map.mapActive.genMinAmp4th = min(bsArray(fourthWindowStartInd:fourthWindowEndInd,:));
handles.data.map.mapActive.genMaxAmp4th = max(bsArray(fourthWindowStartInd:fourthWindowEndInd,:));

% THRESHOLDED MIN
result = arrayEventFinder(dataArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.synNegThresh, 'negative');
events = zeros(numTraces,1);
if ~isempty(result)
    tracesWithEvents = result(:,1);
    timeIndicesOfEvents = result(:,2);
    onsetLatEvents = result(:,2)/sr;
    events(tracesWithEvents) = onsetLatEvents;
end
handles.data.map.mapActive.genLatOfMinOnset = events;

% THRESHOLDED MIN -- spontaneous event analysis of baseline interval
result = arrayEventFinder(dataArray(baselineStartIndex:baselineEndIndex,:), handles.data.map.mapActive.synNegThresh, 'negative');
events = zeros(numTraces,1);
if ~isempty(result)
    tracesWithEvents = result(:,1);
%     timeIndicesOfEvents = result(:,2);
%     onsetLatEvents = result(:,2)/sr;
%     events(tracesWithEvents) = onsetLatEvents;
end
handles.data.map.mapActive.spontEventNum = numel(tracesWithEvents);
totalBaselineTime = numTraces * (handles.data.analysis.baselineEnd - handles.data.analysis.baselineStart);
handles.data.map.mapActive.spontEventRate = numel(tracesWithEvents) / totalBaselineTime;
% NEW WAY -- not fully implemented:
% % % array = dataArray(baselineStartIndex : baselineEndIndex, :);
% % % amplitudeThreshold = handles.data.map.mapActive.synNegThresh;
% % % polarity = 'negative';
% % % durationThreshold = synDuration;
% % % result = arrayEventFinderWithDurThresh(array, amplitudeThreshold, polarity, durationThreshold);
% % % events = zeros(numTraces,1);
% % % if ~isempty(result)
% % %     tracesWithEvents = result.tracesWithEvents;
% % %     numberOfEvents = result.numberOfEvents;
% % %     spontEventNum = sum(numberOfEvents);
% % % end
% % % handles.data.map.mapActive.numberOfEvents = numberOfEvents;
% % % handles.data.map.mapActive.spontEventNum = spontEventNum;
% % % totalBaselineTime = numTraces * (handles.data.analysis.baselineEnd - handles.data.analysis.baselineStart);
% % % handles.data.map.mapActive.spontEventRate = numel(tracesWithEvents) / totalBaselineTime;

% THRESHOLDED MAX
% % val = get(handles.eventPolaritySyn, 'Value') % gs 2006 11 19
% % str = get(handles.eventPolaritySyn, 'String')
% % polarity = str{val}
% % result = arrayEventFinder(dataArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.synPosThresh, polarity);
result = arrayEventFinder(dataArray(stimOnInd:synEndInd,:), handles.data.map.mapActive.synPosThresh, 'positive');
events = zeros(numTraces,1);
if ~isempty(result)
    tracesWithEvents = result(:,1);
    timeIndicesOfEvents = result(:,2);
    onsetLatEvents = result(:,2)/sr;
    events(tracesWithEvents) = onsetLatEvents;
end
handles.data.map.mapActive.genLatOfMaxOnset = events;

% TYPES OF RESPONSES according to clamp mode, response polarity
if handles.data.analysis.currentClamp == 0 %V-clamp
    % assumes neg Vh; inward currents excitatory, outward currents inhibitory

    % first the direct responses: suprathreshold excitatory responses starting within direct window
    handles.data.map.mapActive.edBinary = handles.data.map.mapActive.genLatOfMinOnset > (directWindowStart - stimOn) & ...
        handles.data.map.mapActive.genLatOfMinOnset <= directWindow;
    handles.data.map.mapActive.edAmp = handles.data.map.mapActive.genMinAmp.*handles.data.map.mapActive.edBinary;
    handles.data.map.mapActive.edDelay = handles.data.map.mapActive.genLatOfMin.*handles.data.map.mapActive.edBinary;
    handles.data.map.mapActive.edDelayOnset = handles.data.map.mapActive.genLatOfMinOnset.*handles.data.map.mapActive.edBinary;

    % now the synaptic responses -- suprathreshold excitatory responses starting within synaptic window
    handles.data.map.mapActive.ecBinary = handles.data.map.mapActive.genLatOfMinOnset > synapticWindowStart & ...
        handles.data.map.mapActive.genLatOfMinOnset <= synapticWindow;
    handles.data.map.mapActive.ecAmp = handles.data.map.mapActive.genMinAmp.*handles.data.map.mapActive.ecBinary;
    handles.data.map.mapActive.ecDelay = handles.data.map.mapActive.genLatOfMin.*handles.data.map.mapActive.ecBinary;
    handles.data.map.mapActive.ecDelayOnset = handles.data.map.mapActive.genLatOfMinOnset.*handles.data.map.mapActive.ecBinary;

    handles.data.map.mapActive.icBinary = handles.data.map.mapActive.genLatOfMaxOnset > synapticWindowStart & ...
        handles.data.map.mapActive.genLatOfMaxOnset <= synapticWindow;
    handles.data.map.mapActive.icAmp = handles.data.map.mapActive.genMaxAmp.*handles.data.map.mapActive.icBinary;
    handles.data.map.mapActive.icDelay = handles.data.map.mapActive.genLatOfMax.*handles.data.map.mapActive.icBinary;
    handles.data.map.mapActive.icDelayOnset = handles.data.map.mapActive.genLatOfMaxOnset.*handles.data.map.mapActive.icBinary;

elseif handles.data.analysis.currentClamp == 1 % I-clamp mode(s)
    % assumes typical negative Vr -- positive responses excitatory, negative responses inhibitory
    
    % first the direct responses: suprathreshold excitatory responses starting within direct window
    handles.data.map.mapActive.edBinary = handles.data.map.mapActive.genLatOfMaxOnset > (directWindowStart - stimOn) & ...
        handles.data.map.mapActive.genLatOfMaxOnset <= directWindow;
    handles.data.map.mapActive.edAmp = handles.data.map.mapActive.genMaxAmp.*handles.data.map.mapActive.edBinary;
    handles.data.map.mapActive.edDelay = handles.data.map.mapActive.genLatOfMax.*handles.data.map.mapActive.edBinary;
    handles.data.map.mapActive.edDelayOnset = handles.data.map.mapActive.genLatOfMaxOnset.*handles.data.map.mapActive.edBinary;

    % now the synaptic responses -- suprathreshold excitatory responses starting within synaptic window
    handles.data.map.mapActive.epBinary = handles.data.map.mapActive.genLatOfMaxOnset > synapticWindowStart & ...
        handles.data.map.mapActive.genLatOfMaxOnset <= synapticWindow;
    handles.data.map.mapActive.epAmp = handles.data.map.mapActive.genMaxAmp.*handles.data.map.mapActive.epBinary;
    handles.data.map.mapActive.epDelay = handles.data.map.mapActive.genLatOfMax.*handles.data.map.mapActive.epBinary;
    handles.data.map.mapActive.epDelayOnset = handles.data.map.mapActive.genLatOfMaxOnset.*handles.data.map.mapActive.epBinary;

    handles.data.map.mapActive.ipBinary = handles.data.map.mapActive.genLatOfMinOnset > synapticWindowStart & ...
        handles.data.map.mapActive.genLatOfMinOnset <= synapticWindow;
    handles.data.map.mapActive.ipAmp = handles.data.map.mapActive.genMinAmp.*handles.data.map.mapActive.ipBinary;
    handles.data.map.mapActive.ipDelay = handles.data.map.mapActive.genLatOfMin.*handles.data.map.mapActive.ipBinary;
    handles.data.map.mapActive.ipDelayOnset = handles.data.map.mapActive.genLatOfMinOnset.*handles.data.map.mapActive.ipBinary;

end

% ASSIGN responseType based on binary vectors for ed, ec/ep, ic/ip, nr
if handles.data.analysis.currentClamp == 0 % V-clamp
    bitArray = [num2str(handles.data.map.mapActive.icBinary) ...
            num2str(handles.data.map.mapActive.ecBinary) ...
            num2str(handles.data.map.mapActive.edBinary)];
    codes = bin2dec(bitArray);
    translator = ['NR'; 'ED'; 'EC'; '  '; 'IC'; '  '; 'EC'; '  '];
    handles.data.map.mapActive.responseType = translator(codes+1,:);
elseif handles.data.analysis.currentClamp == 1 % I-clamp
%     bitArray = [num2str(handles.data.map.mapActive.ipBinary) ...
%             num2str(handles.data.map.mapActive.epBinary) ...
%             num2str(handles.data.map.mapActive.apBinary)];
    bitArray = [num2str(handles.data.map.mapActive.ipBinary) ...
            num2str(handles.data.map.mapActive.epBinary) ...
            num2str(handles.data.map.mapActive.edBinary)];
    codes = bin2dec(bitArray);
    translator = ['NR'; 'AP'; 'EP'; '  '; 'IP'; '  '; 'EP'; '  '];
    handles.data.map.mapActive.responseType = translator(codes+1,:);
end

% STORE map arrays for this map
handles = makeMapDataArray(handles);

% PLOT map stuff for this map
handles = array2Dplot(handles);
